Muscle synergies inherent in simulated hypogravity running reveal flexible but not unconstrained locomotor control

With human space exploration back in the spotlight, recent studies have investigated the neuromuscular adjustments to simulated hypogravity running. They have examined the activity of individual muscles, whereas the central nervous system may rather activate groups of functionally related muscles, known as muscle synergies. To understand how locomotor control adjusts to simulated hypogravity, we examined the temporal (motor primitives) and spatial (motor modules) components of muscle synergies in participants running sequentially at 100%, 60%, and 100% body weight on a treadmill. Our results highlighted the paradoxical nature of simulated hypogravity running: The reduced mechanical constraints allowed for a more flexible locomotor control, which correlated with the degree of spatiotemporal adjustments. Yet, the increased temporal (shortened stance phase) and sensory (deteriorated proprioceptive feedback) constraints required wider motor primitives and a higher contribution of the hamstring muscles during the stance phase. These results are a first step towards improving astronaut training protocols.

; GM, gluteus maximus; VM, vastus medialis; VL, vastus lateralis; RF, rectus femoris; STSM, semitendinosus/ semimembranosus; BF, biceps femoris; GaM, gastrocnemius medialis; GaL, gastrocnemius lateralis; SOL, soleus; TA, tibialis anterior; PL, peroneus longus) were collected from healthy men (n = 38) during 2 runs (RUN1 and RUN2) at preferred speed on a LBPPT.Each run consisted of 3 conditions (INIT, HYPO, RLD).The darker rectangles indicate the 60 running cycles considered at the end of each condition.(b) Muscle synergies were extracted from pre-processed sEMG signals using NMF.They consisted of time-invariant motor modules activated by time-varying motor primitives, whose linear combination approximated the pre-processed sEMG signals.(c) Motor primitives were compared across conditions and runs using linear and nonlinear metrics.The center of activity (CoA) indicated when their main activation occurred in time and the full width at half maximum (FWHM) indicated their duration.The HFD described their local complexity (i.e.roughness within each running cycle) and the HE described their global complexity (i.e.irregularity across running cycles).The explanatory diagrams for the nonlinear metrics are adapted from the work of Santuz et al. 34 .
K-means clustering identified four muscle synergies across all three conditions (Fig. 4).The first synergy was functionally related to the braking phase: Its main activation occurred during the first half of the stance phase, with a major contribution from the hip and knee extensors muscles (GM, VM, VL and RF).The second synergy described the push-off phase: Its main activation occurred during the second half of the stance phase, with a major contribution from the ankle plantarflexor muscles (GaM, GaL, SOL and PL).The third synergy was related to the early flight phase, with a major contribution from the ankle dorsiflexor (TA) and evertor (PL) muscles.The fourth synergy differed greatly between conditions: In INIT and RLD, it was functionally related to the late flight phase (preactivation before impact), whereas in HYPO it presented an additional activation during the push-off phase.In all conditions, it showed a major contribution from the hamstring muscles (BF and STSM).In the following, we will refer to this fourth synergy as "S4".Note that some participants did not show all four muscle synergies mentioned above (Supplementary Table 3).www.nature.com/scientificreports/

Increased hamstring contribution
The motor modules were first compared between INIT and HYPO by examining the muscle contributions (Fig. 4a).During the braking phase, the SOL contribution decreased (−22 ± 27%, ES = 0.5), but the GM contribution increased (+29 ± 18%, ES = 0.5).No difference was found in the muscle contributions during the pushoff phase.However, the S4 motor primitives showed an additional activation during this phase (see Fig. 4c for detailed analysis), highlighting a new contribution of the hamstring muscles (STSM and BF) to this phase.During the early flight phase, the RF and GM contributions decreased (RF: −51 ± 17%, ES = 0.7 and GM: −23 ± 27%, ES = 0.3), but the hamstring contribution increased (STSM: +268 ± 12%, ES = 1.1 and BF: +795 ± 18%, ES = 1.0).The S4 motor modules showed a decreased TA contribution (−67 ± 16%, ES = 0.9).They were then compared by calculating a co-contribution index (CI), which assessed the proportion of contributions from the anterior and posterior lower limb muscles at the three joints considered (hip, knee and ankle) (Fig. 5, see Methods for details).From INIT to HYPO, the CI decreased at the hip and knee joints during the early flight phase (−37 ± 6%, ES = 1.0 and −48 ± 5%, ES = 1.3, respectively), indicating an increased contribution of the hamstrings at the expense of the quadriceps.It also decreased at the ankle joint in the S4 motor modules (−46 ± 6%, ES = 1.1), indicating an increased contribution of the triceps surae at the expense of the TA.

Widening of stance motor primitives
Motor primitives were compared across conditions using linear metrics.Specifically, their center of activity (CoA) and full width at half maximum (FWHM) were calculated.They provided an indication of when their main activation occurred in the running cycle and their duration, respectively (Fig. 4b).From INIT to HYPO, the push-off and early flight motor primitives were delayed in the running cycle (higher CoA; +20 ± 23%, ES = 1.3 and +9 ± 14%, ES = 0.7, respectively).Furthermore, the braking and push-off motor primitives were of longer duration relative to the running cycle (higher FWHM; +17 ± 35%, ES = 0.8 and +13 ± 35%, ES = 0.8, respectively), while the early flight primitive was of shorter duration (lower FWHM; −20 ± 26%, ES = 0.4).
The S4 motor primitives differed greatly between INIT and HYPO, with their additional activation in the push-off phase giving them a bimodal profile in HYPO.As they could not be compared using classical metrics, they were excluded from the CoA and FWHM analysis.Instead, the S4 motor primitives were compared across conditions using one-dimensional statistical parametric mapping 35 (SPM) (Fig. 4c).This allowed the calculation of the t-statistics, denoted SPM{t}, for each time point of the running cycle (from 1 to 200).When SPM{t} exceeded the critical threshold computed using random field theory (t*), the motor primitives were considered to be significantly different.From INIT to HYPO, the SPM analysis mainly highlighted the aforementioned additional activation of the S4 motor primitives during the push-off phase (between time points 52-122).
The frequency at which the motor primitives exceeded half of their maximum of the 60 running cycles considered was computed for each time point of the running cycle (Fig. 6a).This enabled to calculate the frequency      of overlaps (Fig. 6b), which we compared between conditions using SPM (Fig. 6c): An overlap occurs when at least two primitives simultaneously exceed half of their maximum.Consistent with the aforementioned widening of the braking and push-off motor primitives, the heat maps showed that both exceeded half of their maximum over a larger range of the running cycle in HYPO than in INIT and RLD.However, these adjustments were not systematic over the 60 running cycles considered, as shown by the lighter shades of grey and red on the heat maps.This resulted in a decreased frequency of overlaps between the braking and push-off motor primitives (between time points 22-34).On the other hand, the heat maps highlighted the additional activation of the S4 motor primitives during the push-off phase in HYPO, leading to an increased frequency of overlaps between the push-off and S4 motor primitives (between time points 48-81).www.nature.com/scientificreports/

Correlations between adjustments in the spatiotemporal parameters of running and motor primitives
Finally, the adjustments in the spatiotemporal parameters of running (stance time, flight time, stride frequency) were correlated with those in motor primitives (FWHM, COA, HFD, HE) (Fig. 8).When focusing on linear are color-coded: from white (the primitives never exceeded half of their maximum) to the darkest color (the primitives exceeded half of their maximum in all running cycles considered).Missing motor primitives are reported as fully white rows (participants did not systematically present the four muscle synergies, see Supplementary Table 3).metrics, correlation analysis showed that the participants with the greatest decrease in stance time had the greatest increase in the FWHM of the push-off primitives (ρ = −0.41,p < 0.05).Furthermore, those with the greatest increase in flight time had the greatest increase in the FWHM of the braking primitives and in the CoA of the push-off primitives (ρ = 0.40 and ρ = 0.35, p < 0.05 for both) (Fig. 8a).When looking at the non-linear metrics derived from fractal theory, correlational analysis showed that the participants with the smallest decrease in stance time had the greatest increase in the HFD of the push-off primitives (ρ = 0.51, p < 0.05).Finally, those with the greatest decrease in stride frequency had the greatest increase in the HFD of the braking and S4 primitives (ρ = −0.33 and ρ = −0.44,p < 0.05 for both), and in the HE of the S4 primitives (ρ = −0.52,p < 0.05) (Fig. 8b).

Discussion
The present study provided a better understanding of locomotor control during simulated hypogravity running.As a preamble, simulated hypogravity running resulted in the adoption of a bouncing running pattern, characterized by a slowing of stride frequency.Although the gravity simulated in this study (0.6 g) is closer to that of Mars than that of the Moon, these biomechanical adjustments are largely consistent with videos from the Apollo 11-17 missions, which showed that the astronauts neither walked nor ran, but rather adopted hopping, skipping 36 or loping 37 gaits.This suggests a specific reorganization of locomotor control, which we investigated by analyzing the muscle synergies in the light of fractal theory.The extracted muscle synergies described the four usual phases of the running cycle, both in normogravity and in simulated hypogravity.In this sense, it has been reported that the postural and walking muscle synergies were consistent up to 0.25 g and 0.16 g, respectively 38,39 .This suggests that the training protocols that astronauts undergo in normogravity (on Earth) may be transferable to hypogravity (e.g., on the Moon or Mars) and vice versa.However, it should be noted that three muscle synergies would be sufficient to account for the sEMG signals when walking in very low gravity (0.07 g), compared to four at higher gravities (from 1 to 0.16 g) 39 .Although the current study only examined the running muscle synergies at 0.6 g, which is far from the gravity of Mars and the Moon, it highlighted some adjustments in the motor modules and motor primitives that should be considered in anticipation for future space missions.
Simulated hypogravity running led to a redistribution of muscle contributions within the motor modules.During the braking phase, the reduced SOL contribution is consistent with its antigravity function 40 .In the S4 motor modules, the decreased TA contribution is attributed to the shift toward a forefoot strike pattern 41,42 .Importantly, the push-off phase, showed no change in the contribution of the ankle plantarflexors (SOL, GaM, GaL, PL), but an increased contribution of the hamstring muscles (BF and STSM), as evidenced by the additional activation of the S4 motor primitive.This was previously reported during simulated hypogravity walking 43 and confirms our recent findings of increased hamstring activity during simulated hypogravity running, especially in the braking (BF: + 44 ± 18%) and push-off (BF: + 49 ± 12%, STSM: + 123 ± 14%) phases 11 .This may have compensated for the disadvantageous force output of the gastrocnemii in a forefoot versus rearfoot strike pattern 44 .www.nature.com/scientificreports/ In addition, the hamstring muscle group is known to be involved in increasing walking and running speed 45,46 .Its increased contribution may thus have contributed to accelerating the knee joint flexion of the swing leg, in an attempt to counteract the slowing of stride frequency, characteristic of simulated hypogravity running.This is supported by its increased contribution in the early flight phase, which came at the expense of the RF contribution due to reciprocal inhibition 47 .However, the adjustments in hamstring contribution were not correlated with those in stride frequency.Further studies should therefore investigate the functional role of this muscle group during simulated hypogravity running in order to provide recommendations for pre-flight training of astronauts.
As things stand, specific strengthening of this muscle group may be recommended.Simulated hypogravity running resulted in wider stance motor primitives (higher FWHM of the braking and push-off primitives).Furthermore, the push-off primitives were shifted later in the running cycle (higher CoA).Such adjustments of the stance primitives has previously been observed in genetically modified mice lacking muscle spindles and Golgi tendon organs when compared to wild type 34,48 .These results could therefore be attributed to the deteriorated proprioceptive feedback caused by simulated hypogravity 49,50 .Yet, the opposite adjustments found for early flight motor primitives, which were narrower, rather suggests the influence of temporal constraints: The stance phase was shortened during simulated hypogravity running, with the push-off phase being more affected than the braking phase 11 , while the flight phase was prolonged.This is supported by the correlation analysis, which showed that the participants with the greatest decrease in stance time were those with the greatest widening of the push-off motor primitives.Thus, the lack of time to organize locomotor control during the stance phase may have been compensated for by increasing the overlap between the temporally adjacent braking and push-off motor primitives.Typically reported in the presence of external perturbations [26][27][28][29][30] , such adjustment was not systematically observed across running cycles, as indicated by the decreased frequency  3).The blue shaded areas correspond to the 95% confidence interval.Non-significant correlations are shown in Supplementary Table 6 www.nature.com/scientificreports/ of overlaps.This is a first indication that simulated hypogravity running may constrain locomotor control, while allowing its flexibility.The nonlinear analysis of the motor primitives supported this assumption of a higher locomotor control flexibility.On the one hand, three of the four primitives (those of braking, early flight, and S4) were more locally complex (higher HFD).This result may seem surprising given the literature on externally perturbed locomotion, which reports a decrease in the HFD of motor primitives 27,28 .However, it was expected since simulated hypogravity running is certainly perturbed, but also less mechanically contrained.Only the push-off primitives were not more locally complex.This is again attributed to the aforementioned high temporal constraints on push-off phase 11 .Indeed, the correlation analysis showed that participants with the greatest decrease in stance time were those showing the lowest increase in the HFD of the push-off motor primitives.On the other hand, all motor primitives were more globally complex (higher HE) during simulated hypogravity running.Since the HE increased from 0 (sinusoidal time series) to 0.5 (random time series), it follows that the motor primitives were more irregular across running cycles in simulated hypogravity than in normogravity.This suggests that participants exerted less stride-to-stride control 51 .The unintended horizontal support provided by the LBPPT chamber 52 , in addition to the expected vertical support, could explain this reduced control effort.It would therefore be relevant to study its influence on muscle synergies.In any case, this is in line with the 'minimum intervention principle' used by humans to regulate locomotion: Deviations are corrected only when they interfere with task performance 53 .
Is this increase in local and global complexity of motor primitives beneficial for simulated hypogravity running?Correlation analysis showed that participants with the greatest increase in the HFD of the braking and S4 motor primitives, and in the HE of the S4 motor primitives had the greatest decrease in stride frequency.All correlations were moderate to large, indicating a non-negligible relationship between these nonlinear metrics and stride frequency.Because it is systematically modified, stride frequency appears to be the spatiotemporal parameter that best describes the degree of adjustment to simulated hypogravity running.In this sense, it has been proposed that in mature and healthy biological systems, an optimal amount of variability reflects the flexibility of locomotor control to perturbations 32,33 .Less than optimal variability makes the system overly rigid, while greater than optimal variability makes the system overly noisy and unstable.Thus, the higher local and global complexity of the motor primitives observed in this study would reflect a more flexible locomotor control, allowing participants to optimally adjust to simulated hypogravity running.This finding may pave the way for the development of a feedback tool to improve astronauts' ability to adjust to hypogravity.
Overall, our results revealed a specific reorganization of locomotor control in the face of the paradoxical nature of simulated hypogravity running, which combines reduced mechanical constraints with increased temporal and sensory constraints.With human space exploration back in the spotlight, they deserve to be verified in simulated gravities closer to those of Mars or the Moon.They also need be confirmed in actual hypogravity: Subtle but systematic differences in the running pattern are known, suggesting potential differences in locomotor control 54 .In any case, the current results provide a first insight into locomotor control during simulated hypogravity running.They are a first step towards improving pre-flight training, per-flight countermeasures and post-flight rehabilitation of astronauts.

Participants
Thirty-eight recreational runners (age: 19 ± 1 years old [mean ± sd]; height: 177 ± 6 cm; mass: 69.3 ± 8.1 kg) volunteered for this study.We ensured that they had not suffered any musculoskeletal injury for more than 1 year, and had no history of any neurological impairment.All procedures were approved by the National Ethics Committee for Research in Sports Sciences (reference number: CERSTAPS IRB00012476-2021-31-03-96).All methods were carried out in accordance with the Declaration of Helsinki, and written informed consent was obtained from all participants prior to the experiment.

Experimental design
Participants were asked to run on a Lower Body Positive Pressure Treadmill (LBPPT) (VIA400X AlterG®, Fremont, CA, USA) at their preferred speed (3.0 ± 0.2 m s -1 ).They were instructed to remain in the center of its chamber to minimize involuntary horizontal support 52 .The preferred speed was self-selected while running at 100% body weight (1 g) during a familiarization session performed the week pior to the experimental protocol.After a 3-min warm-up to gradually reach the preferred speed, the experimental protocol consisted of two 9-min runs separated by a 4-min recovery period.Each run included 3 consecutive 3-min running conditions: The initial (INIT) and reloaded (RLD) conditions were run at 100% body weight, while the intermediate simulated hypogravity condition (HYPO) was run at 60% body weight (0.6 g).The transitions phases between the conditions were progressive and lasted 12 ± 2 s.

Spatiotemporal parameters of running
Touchdown and toe-off were identified from the normal component of the ground reaction force using a threshold of 50 N.Stance time, flight time and stride frequency were calculated.The foot strike pattern (rearfoot or midfoot/forefoot) was obtained from sensors located in the anterior and posterior regions of the instrumented insoles.

EMG pre-processing
The raw sEMG signals were first band-pass filtered (cutoff at 10 and 500 Hz).They were then high-pass filtered (cutoff at 50 Hz), full-wave rectified, and low-pass filtered (cutoff at 20 Hz) using a 4th order IIR Butterworth zero-phase filter 24 .After subtracting the minimum, the sEMG signals were normalized to the maximum activity recorded for every muscle in each running condition 56 .Each running cycle was then time-normalized to 200 points, assigning 100 points to the stance phase and 100 points to the flight phase.

Muscle synergy extraction
Muscle synergies were extracted from the pre-processed sEMG signals using a custom R script (v3.6.3,R Foundation for Statistical Computing, Vienna, Austria).The m = 11 muscles listed above were included in the analysis (TA, GaM, GaL, SOL, PL, VM, VL, RF, BF, STSM, GM).For each condition and participant, the m = 11 timedependent muscle activity vectors were grouped into a matrix V with dimensions m = 11 rows and n = 1200 columns (200 points × 60 running cycles).The matrix V was factorized using the classical Gaussian Non-negative Matrix Factorization (NMF) algorithm 22,57 such that V ≈ V R = WH.The new matrix V R , reconstructed by multiplying the matrices W and H, approximates the initial matrix V.The motor modules matrix W contained the time-invariant muscle weights that describe the relative contribution of each muscle within a given synergy.It had dimensions m × r, where "r" is the minimum number of synergies necessary to satisfactorily reconstruct the original matrix V.The motor primitive matrix H, with dimensions s × n, contained the time-dependent activation of each synergy.The update rules for W and H are presented in Eqs. ( where "T" is the transposed matrix and ⊙ the element-by-element multiplication.The quality of the reconstruction was assessed by calculating the coefficient of determination (R 2 ) between the original and reconstructed matrices (V and V R , respectively), as shown in Eqs. ( 3) and (4).
where SSE is the sum squared error, SST the sum squared total and V the average of the matrix V.
Iterations were stopped when the change in the calculated R 2 was less than 0.01% over the last 20 iterations 24 .This operation was first completed by setting the number of synergies to 1.It was then repeated by increasing the number of synergies, up to a maximum of 8 synergies.The number of 8 was chosen to be equal to 75% of the number of monitored muscles, as extracting a number of synergies equal to the number of monitored muscles would not reduce the dimensionality of the data.For each number of synergies, the factorization was repeated 5 times, each time generating new randomized initial matrices W and H, to avoid local minima 58 .The solution with the highest R 2 was then selected.To select the minimum number of synergies required to reconstruct the original sEMG signals, the curve of R 2 values vs. number of synergies was fitted using a simple linear regression model.The mean squared error between the curve and the linear regression was calculated 59 .The lowest abscissa point was then removed and the error between this new curve and its new linear regression was calculated.This was repeated until there were 2 points left on the curve or the mean squared error was less than 10 -4 .This was done to find the most linear part of the curve of R 2 values vs. number of synergies, assuming that in this section the reconstruction quality could not be significantly increased by adding more synergies to the model.

Co-contribution index
For each lower limb joint (hip, knee, ankle), we calculated a co-contribution index (CI) from the motor modules.To do this, we computed the mean contributions of the anterior ( ant ) and posterior ( post ) mobilizing muscles 60 .For the hip, the anterior mobilizing muscle was the RF anterior and the posterior was the GM.For the knee, the anterior mobilizing muscles were RF, VL and VM and the posterior were the STSM and BF.For the ankle, (1) www.nature.com/scientificreports/ the anterior mobilizing muscle was the TA and the posteriors were the GaM, GaL, SOL and PL.Equation ( 5) defines the CI.
According to this definition, (i) CI = 0 when only the posterior muscles contribute to the considered joint; (ii) CI = 1 when only the anterior muscles contribute and (iii) CI = 0.5 if anterior and posterior muscles contribute equally.

Center of activity of motor primitives
We compared the averaged motor primitives (over the 60 considered running cycles) by evaluating first their center of activity (CoA).The CoA is an indication of when their main activation happens in time.It was calculated as the angle of the vector (in polar coordinates) oriented towards the center of mass of that circular distribution 61 .The polar direction represented the running cycle's phase, 0° corresponding to the touchdown and 360° to the next touchdown.Equations ( 6), ( 7) and ( 8) define the CoA.
where "p" is the number of points of the running cycle (p = 200) and H is the motor primitive averaged over the 60 running cycles.

Width of motor primitives
We also compared the averaged motor primitives based on their full width at half maximum (FWHM), which describe their duration.It was calculated, after subtracting their minimum, as the number of points exceeding half of their maximum 25,26 .For each point of the running cycle (from 1 to 200), we also calculated the frequency at which motor primitives exceeded half of their maximum over the 60 running cycles, and the frequency at which they overlapped.Overlap occurred when at least two primitives exceeded half of their maximum at the same time 30,62 .

Local complexity of motor primitives
To assess the local complexity of motor primitives (i.e.their "roughness" within each running cycle), we calculated their Higuchi's fractal dimension 63 (HFD).For each motor primitive H(t) = [H(1) H(2) … H(n)], we constructed k new times series as described in Eq. ( 9), where k ∈ ⟦2 ; k max ⟧.
where t 0 ∈ ⟦1; k⟧ is the first sample at initial time and ⌊ ⌋ the integer part.Then, we calculated the non-Euclidean Higuchi length L t 0 k of each new time series H t 0 k as defined in Eq. (10).
The length of the motor primitive L k was defined as the average of the k sets lengths as shown in Eq. ( 11).
If L k ∝ k −HFD , then the curve is fractal with dimension HFD and this should lead the plot of log(L k ) vs. log(1/k) to fall on a straight line with slope HFD.The subsampling rate k was determined by repeating the pre- vious operation for k = 1, 2, 3, … up to a maximum of k = 300, and then selecting the maximum value of k, k max , for which the curve of log(L k ) vs. log(1/k) was linear 31 .Our data led us to choose k max = 10.HFD varies between 1 and 2, with increasing values corresponding to a more locally complex (or rough) signal, and the value 1.5 indicating a random Gaussian noise 63,64 .

Global complexity of motor primitives
To assess the global complexity of motor primitives (i.e.their irregularity across running cycles), we calculated their Hurst exponent (HE) following the rescaled range approach 65,66 .For each motor primitive H(t) = [H(1) H(2) … H(n)], we constructed a new time series of length q (q ∈ ⟦1 ; n⟧), as described in Eq. ( 12).
H q was then centred by substracting its mean H q and a series of cumulative values H t q was established, as described in Eq. 13.
We determined the range R q of the time series H t q (i.e. the difference between the maximum and minimum values) and its standard deviation S q .We computed its rescaled range (R q /S q ) .The value of q was determined by repeating the previous operations for q = n, n/2, n/4, n/8, … until a minimum of q = 2, and then selecting the value of q, q inflex , where an inflexion point occurred 31 .We finally calculated HE as the slope of the curve of log q vs. log R q /S q between q = 2 and q = q inflex .Our data led us to choose q inflex = 200.HE can vary between 0 and 1.For 0.5 < HE < 1, the time series (i.e.motor primitive) is persistent, meaning that high values will probably be followed by other high values (and vice versa).For 0 < HE < 0.5, the time series is anti-persistent meaning that high values will probably be followed by low values (with frequent alternation between the two).HE = 0.5 indicates a completely random series without any persistence 66 .

Functional classification of muscle synergies
All the above calculations were conducted only for the motor primitives relative to fundamental synergies.The recognition of fundamental synergies was carried out by k-means clustering to reduce possible operator bias in the classification.We first clustered the average motor primitives (over the 60 running cycles), for each condition separately.This classification operation was performed for a number of clusters ranging from 1 to 11 using the Hartigan Wong algorithm 67 .For each number of clusters, the classification was repeated 20 times with different initializations of the centroids.The solution that minimized the intra-class variance was retained.The number of clusters was determined by fitting the curve of intra-class variance vs. number of clusters with a simple linear regression model, and processing in the same way as for the determination of the number of synergies (until mean squared error was less than 10 -4 ).The motor modules were then classified by imposing the same number of clusters as previously obtained for the motor primitives.The average FWHM and CoA of the motor primitives were then summed and normalized by the number of points (i.e.200).This value was used as a score to compare the k-means clustering of motor primitives and motor modules.Consistent clusters identified fundamental synergies, while discordant clusters identified combined synergies (linear combination of other simpler synergies).Due to the lack of consensus in the literature on how to interpret combined synergies, they were excluded from the analysis.

Statistics
For each variable, a linear mixed model (R package LmerTest 68 ) was performed using restricted maximum likelihood estimation.Joint, Muscle, Synergy, Condition and Run were considered as fixed effects when appropriate, and preferred running speed as covariate.The intercepts for the participants and the slope per condition depending on the foot strike pattern were chosen as random effects.The significance of the random effects was tested.The best model (i.e. the number of fixed effects and interactions) was selected by likelihood ratio tests of model comparisons using a backward selection method.An Analysis of Variance (ANOVA) (degrees of freedom estimated with the Satterthwaite formula) was then performed on the selected model.This was followed by multiple comparisons with Tukey (for main effects) and Holm (for interaction effects) adjustment.If the normality of the residuals was violated (Shapiro-Wilk test), a permutation ANOVA of the linear mixed model was performed (R package lmPerm, number of permutations: 10,000).The effect size (ES) was calculated using Cohen's d coefficient 69 and assessed using the following thresholds: 0.2 to < 0.6, 0.6 to < 1.2 and 1.2 to < 2.0 for small, moderate and large effects, respectively 70 .In cases where motor primitives could not be compared using conventional linear metrics (CoA and FWHM) because they exhibited a bimodal profile, we used a two-way (Condition, Run) repeated measures ANOVA based on the one-dimensional Statistical Parametric Mapping analysis (SPM) (Python package spm1d) 35 .This was followed by multiple comparisons with Bonferroni adjustment.Finally, Spearman's correlation tests were used to quantify the correlations between the adjustments in motor primitives (in CoA, FWHM, HFD and HE) and those in the spatiotemporal parameters of the running pattern (in stance time, flight time and stride frequency).The magnitude of Spearman correlation coefficients was interpreted using the following thresholds: 0.1 to < 0.3, 0.3 to < 0.5 and > 0.5 for small, moderate and large, respectively 70 .
All the significance levels were set to α = 0.05 and the statistical analyses were performed using custom R and Python scripts (v3.9.7, 2021, Python Software Foundation, Wilmington, Delaware, USA).All results are presented in the text as estimated mean ± standard error.

Data availability
The raw data and the codes developed for their analysis are available to any researcher upon request to the corresponding author: Camille Fazzari (camille.fazzari@univ-amu.fr).( 12) H q = [H(1) H(2) . . .H q ] (13) H t q = t u=1 H q (u) − H q

Figure 1 .
Figure 1.Experimental design and muscle synergy analyses.(a) Normal ground reaction force and sEMG activity (11 muscles; GM, gluteus maximus; VM, vastus medialis; VL, vastus lateralis; RF, rectus femoris; STSM, semitendinosus/ semimembranosus; BF, biceps femoris; GaM, gastrocnemius medialis; GaL, gastrocnemius lateralis; SOL, soleus; TA, tibialis anterior; PL, peroneus longus) were collected from healthy men (n = 38) during 2 runs (RUN1 and RUN2) at preferred speed on a LBPPT.Each run consisted of 3 conditions (INIT, HYPO, RLD).The darker rectangles indicate the 60 running cycles considered at the end of each condition.(b) Muscle synergies were extracted from pre-processed sEMG signals using NMF.They consisted of time-invariant motor modules activated by time-varying motor primitives, whose linear combination approximated the pre-processed sEMG signals.(c) Motor primitives were compared across conditions and runs using linear and nonlinear metrics.The center of activity (CoA) indicated when their main activation occurred in time and the full width at half maximum (FWHM) indicated their duration.The HFD described their local complexity (i.e.roughness within each running cycle) and the HE described their global complexity (i.e.irregularity across running cycles).The explanatory diagrams for the nonlinear metrics are adapted from the work of Santuz et al.34 .

Figure 3 .
Figure 3. Pre-processed sEMG signals.The x-axis distinguishes the stance phase from the flight phase, normalized to the same number of points.The y-axis shows the normalized sEMG activity.Data are averaged over all participants (n = 38), with shaded standard deviation.

Figure 4 .
Figure 4. Comparison of muscle synergies across conditions.(a) Motor modules.The y-axis shows the normalized muscle contributions, with individual data displayed.A (Condition x Synergy x Muscle) interaction was found (permutation ANOVA, F (60, 8491) = 4.6, p < 0.001).(b) Motor primitives.The x-axis distinguishes the stance phase from the flight phase, normalized to the same number of points.The y-axis shows the normalized amplitude with shaded standard deviation.(Condition x Synergy) interactions were found for the FWHM (permutation ANOVA, F (4, 576) = 7.6, p < 0.001) and the CoA (ANOVA, F (4, 576) = 5.1, p < 0.001).c S4 could not be compared using classical metrics (COA and FWHM), but SPM revealed a main effect of the condition (SPM ANOVA, F* = 8.3).The blue shaded supra-threshold clusters (t* = 4.02) correspond to the time periods in which the S4 motor primitives differed.Motor modules and primitives are averaged across participants presenting the synergy of interest (see Supplementary Table3).Significant adjustments determined by post-hoc analysis (multiple comparisons with Holm's adjustment for panels a-b and with Bonferroni's adjustment for panel c) are shown for HYPO vs. INIT (*) and RLD vs. HYPO ($) (see Supplementary Table4-5 for details).The number of symbols indicates the statistical level (one for p < 0.05, two for p < 0.01, three for p < 0.001).

3
Figure 4. Comparison of muscle synergies across conditions.(a) Motor modules.The y-axis shows the normalized muscle contributions, with individual data displayed.A (Condition x Synergy x Muscle) interaction was found (permutation ANOVA, F (60, 8491) = 4.6, p < 0.001).(b) Motor primitives.The x-axis distinguishes the stance phase from the flight phase, normalized to the same number of points.The y-axis shows the normalized amplitude with shaded standard deviation.(Condition x Synergy) interactions were found for the FWHM (permutation ANOVA, F (4, 576) = 7.6, p < 0.001) and the CoA (ANOVA, F (4, 576) = 5.1, p < 0.001).c S4 could not be compared using classical metrics (COA and FWHM), but SPM revealed a main effect of the condition (SPM ANOVA, F* = 8.3).The blue shaded supra-threshold clusters (t* = 4.02) correspond to the time periods in which the S4 motor primitives differed.Motor modules and primitives are averaged across participants presenting the synergy of interest (see Supplementary Table3).Significant adjustments determined by post-hoc analysis (multiple comparisons with Holm's adjustment for panels a-b and with Bonferroni's adjustment for panel c) are shown for HYPO vs. INIT (*) and RLD vs. HYPO ($) (see Supplementary Table4-5 for details).The number of symbols indicates the statistical level (one for p < 0.05, two for p < 0.01, three for p < 0.001).
Figure 4. Comparison of muscle synergies across conditions.(a) Motor modules.The y-axis shows the normalized muscle contributions, with individual data displayed.A (Condition x Synergy x Muscle) interaction was found (permutation ANOVA, F (60, 8491) = 4.6, p < 0.001).(b) Motor primitives.The x-axis distinguishes the stance phase from the flight phase, normalized to the same number of points.The y-axis shows the normalized amplitude with shaded standard deviation.(Condition x Synergy) interactions were found for the FWHM (permutation ANOVA, F (4, 576) = 7.6, p < 0.001) and the CoA (ANOVA, F (4, 576) = 5.1, p < 0.001).c S4 could not be compared using classical metrics (COA and FWHM), but SPM revealed a main effect of the condition (SPM ANOVA, F* = 8.3).The blue shaded supra-threshold clusters (t* = 4.02) correspond to the time periods in which the S4 motor primitives differed.Motor modules and primitives are averaged across participants presenting the synergy of interest (see Supplementary Table3).Significant adjustments determined by post-hoc analysis (multiple comparisons with Holm's adjustment for panels a-b and with Bonferroni's adjustment for panel c) are shown for HYPO vs. INIT (*) and RLD vs. HYPO ($) (see Supplementary Table4-5 for details).The number of symbols indicates the statistical level (one for p < 0.05, two for p < 0.01, three for p < 0.001).

4
Figure 4. Comparison of muscle synergies across conditions.(a) Motor modules.The y-axis shows the normalized muscle contributions, with individual data displayed.A (Condition x Synergy x Muscle) interaction was found (permutation ANOVA, F (60, 8491) = 4.6, p < 0.001).(b) Motor primitives.The x-axis distinguishes the stance phase from the flight phase, normalized to the same number of points.The y-axis shows the normalized amplitude with shaded standard deviation.(Condition x Synergy) interactions were found for the FWHM (permutation ANOVA, F (4, 576) = 7.6, p < 0.001) and the CoA (ANOVA, F (4, 576) = 5.1, p < 0.001).c S4 could not be compared using classical metrics (COA and FWHM), but SPM revealed a main effect of the condition (SPM ANOVA, F* = 8.3).The blue shaded supra-threshold clusters (t* = 4.02) correspond to the time periods in which the S4 motor primitives differed.Motor modules and primitives are averaged across participants presenting the synergy of interest (see Supplementary Table3).Significant adjustments determined by post-hoc analysis (multiple comparisons with Holm's adjustment for panels a-b and with Bonferroni's adjustment for panel c) are shown for HYPO vs. INIT (*) and RLD vs. HYPO ($) (see Supplementary Table4-5 for details).The number of symbols indicates the statistical level (one for p < 0.05, two for p < 0.01, three for p < 0.001).
Figure 4. Comparison of muscle synergies across conditions.(a) Motor modules.The y-axis shows the normalized muscle contributions, with individual data displayed.A (Condition x Synergy x Muscle) interaction was found (permutation ANOVA, F (60, 8491) = 4.6, p < 0.001).(b) Motor primitives.The x-axis distinguishes the stance phase from the flight phase, normalized to the same number of points.The y-axis shows the normalized amplitude with shaded standard deviation.(Condition x Synergy) interactions were found for the FWHM (permutation ANOVA, F (4, 576) = 7.6, p < 0.001) and the CoA (ANOVA, F (4, 576) = 5.1, p < 0.001).c S4 could not be compared using classical metrics (COA and FWHM), but SPM revealed a main effect of the condition (SPM ANOVA, F* = 8.3).The blue shaded supra-threshold clusters (t* = 4.02) correspond to the time periods in which the S4 motor primitives differed.Motor modules and primitives are averaged across participants presenting the synergy of interest (see Supplementary Table3).Significant adjustments determined by post-hoc analysis (multiple comparisons with Holm's adjustment for panels a-b and with Bonferroni's adjustment for panel c) are shown for HYPO vs. INIT (*) and RLD vs. HYPO ($) (see Supplementary Table4-5 for details).The number of symbols indicates the statistical level (one for p < 0.05, two for p < 0.01, three for p < 0.001).

Figure 5 .
Figure 5. Co-contribution from anterior and posterior lower limb muscles.For each joint, the anterior (purple) and posterior (green) muscles considered are shown in the header.Arrows pointing to the right indicate a greater contribution from anterior muscles while arrows pointing to the left indicate a greater contribution from posterior muscles.Arrows pointing upwards indicate a perfect balance between the contribution of the anterior and posterior muscles.An interaction (Condition x Synergy x Joint) was found (permutation ANOVA, F(12, 2313) = 11, p < 0.001).The CI is averaged across participants presenting the synergy of interest (see Supplementary Table3).The standard deviation is shown with curved lines.The significant adjustments determined by post-hoc analysis (multiple comparisons with Holm's adjustment) are shown for HYPO vs. INIT (*) and RLD vs. HYPO ($) (see Supplementary Table5for details).The number of symbols indicates the statistical level (three for p < 0.001).

3
Figure 5. Co-contribution from anterior and posterior lower limb muscles.For each joint, the anterior (purple) and posterior (green) muscles considered are shown in the header.Arrows pointing to the right indicate a greater contribution from anterior muscles while arrows pointing to the left indicate a greater contribution from posterior muscles.Arrows pointing upwards indicate a perfect balance between the contribution of the anterior and posterior muscles.An interaction (Condition x Synergy x Joint) was found (permutation ANOVA, F(12, 2313) = 11, p < 0.001).The CI is averaged across participants presenting the synergy of interest (see Supplementary Table3).The standard deviation is shown with curved lines.The significant adjustments determined by post-hoc analysis (multiple comparisons with Holm's adjustment) are shown for HYPO vs. INIT (*) and RLD vs. HYPO ($) (see Supplementary Table5for details).The number of symbols indicates the statistical level (three for p < 0.001).
Figure 5. Co-contribution from anterior and posterior lower limb muscles.For each joint, the anterior (purple) and posterior (green) muscles considered are shown in the header.Arrows pointing to the right indicate a greater contribution from anterior muscles while arrows pointing to the left indicate a greater contribution from posterior muscles.Arrows pointing upwards indicate a perfect balance between the contribution of the anterior and posterior muscles.An interaction (Condition x Synergy x Joint) was found (permutation ANOVA, F(12, 2313) = 11, p < 0.001).The CI is averaged across participants presenting the synergy of interest (see Supplementary Table3).The standard deviation is shown with curved lines.The significant adjustments determined by post-hoc analysis (multiple comparisons with Holm's adjustment) are shown for HYPO vs. INIT (*) and RLD vs. HYPO ($) (see Supplementary Table5for details).The number of symbols indicates the statistical level (three for p < 0.001).

Figure 6 .
Figure 6.Overlaps of motor primitives.(a)Frequency at which the motor primitives exceeded half of their maximum.Each row of the heatmaps corresponds to a participant, each column to a time point.The results are color-coded: from white (the primitives never exceeded half of their maximum) to the darkest color (the primitives exceeded half of their maximum in all running cycles considered).Missing motor primitives are reported as fully white rows (participants did not systematically present the four muscle synergies, see Supplementary Table3).(b) Frequency of overlaps between motor primitives.The x-axis distinguishes the stance phase from the flight phase, normalized to the same number of points.(c) The frequency of overlaps was compared using SPM, which detected a main effect of the condition (SPM ANOVA, F* = 9.4).The blue shaded supra-threshold clusters (t* = 4.2) correspond to the time periods in which the frequency of overlaps differed.The significant adjustments determined by post-hoc analysis (multiple comparisons with Bonferroni's adjustment) are shown for HYPO vs. INIT (*) and RLD vs. HYPO ($).The number of symbols indicates the statistical level (three forp < 0.001).
Figure 6.Overlaps of motor primitives.(a)Frequency at which the motor primitives exceeded half of their maximum.Each row of the heatmaps corresponds to a participant, each column to a time point.The results are color-coded: from white (the primitives never exceeded half of their maximum) to the darkest color (the primitives exceeded half of their maximum in all running cycles considered).Missing motor primitives are reported as fully white rows (participants did not systematically present the four muscle synergies, see Supplementary Table3).(b) Frequency of overlaps between motor primitives.The x-axis distinguishes the stance phase from the flight phase, normalized to the same number of points.(c) The frequency of overlaps was compared using SPM, which detected a main effect of the condition (SPM ANOVA, F* = 9.4).The blue shaded supra-threshold clusters (t* = 4.2) correspond to the time periods in which the frequency of overlaps differed.The significant adjustments determined by post-hoc analysis (multiple comparisons with Bonferroni's adjustment) are shown for HYPO vs. INIT (*) and RLD vs. HYPO ($).The number of symbols indicates the statistical level (three forp < 0.001).

Figure 7 .
Figure 7. Local and global complexity of motor primitives.(Condition x Synergy) interactions were found for (a) the HFD (ANOVA, F(6,747) = 14, p < 0.001) and (b) the HE (ANOVA, F(6,7450) = 3.4, p < 0.01) of motor primitives.Data are averaged across participants presenting the synergy of interest (see Supplementary Table3), with individual data displayed.Significant adjustments determined by post-hoc analysis (multiple comparisons with Holm's adjustment) are shown for HYPO vs. INIT (*), RLD vs. HYPO ($) (see Supplementary Table5for details).The number of symbols indicates the statistical level (three for p < 0.001).

3
Figure 7. Local and global complexity of motor primitives.(Condition x Synergy) interactions were found for (a) the HFD (ANOVA, F(6,747) = 14, p < 0.001) and (b) the HE (ANOVA, F(6,7450) = 3.4, p < 0.01) of motor primitives.Data are averaged across participants presenting the synergy of interest (see Supplementary Table3), with individual data displayed.Significant adjustments determined by post-hoc analysis (multiple comparisons with Holm's adjustment) are shown for HYPO vs. INIT (*), RLD vs. HYPO ($) (see Supplementary Table5for details).The number of symbols indicates the statistical level (three for p < 0.001).
Figure 7. Local and global complexity of motor primitives.(Condition x Synergy) interactions were found for (a) the HFD (ANOVA, F(6,747) = 14, p < 0.001) and (b) the HE (ANOVA, F(6,7450) = 3.4, p < 0.01) of motor primitives.Data are averaged across participants presenting the synergy of interest (see Supplementary Table3), with individual data displayed.Significant adjustments determined by post-hoc analysis (multiple comparisons with Holm's adjustment) are shown for HYPO vs. INIT (*), RLD vs. HYPO ($) (see Supplementary Table5for details).The number of symbols indicates the statistical level (three for p < 0.001).

5
Figure 7. Local and global complexity of motor primitives.(Condition x Synergy) interactions were found for (a) the HFD (ANOVA, F(6,747) = 14, p < 0.001) and (b) the HE (ANOVA, F(6,7450) = 3.4, p < 0.01) of motor primitives.Data are averaged across participants presenting the synergy of interest (see Supplementary Table3), with individual data displayed.Significant adjustments determined by post-hoc analysis (multiple comparisons with Holm's adjustment) are shown for HYPO vs. INIT (*), RLD vs. HYPO ($) (see Supplementary Table5for details).The number of symbols indicates the statistical level (three for p < 0.001).
Figure 7. Local and global complexity of motor primitives.(Condition x Synergy) interactions were found for (a) the HFD (ANOVA, F(6,747) = 14, p < 0.001) and (b) the HE (ANOVA, F(6,7450) = 3.4, p < 0.01) of motor primitives.Data are averaged across participants presenting the synergy of interest (see Supplementary Table3), with individual data displayed.Significant adjustments determined by post-hoc analysis (multiple comparisons with Holm's adjustment) are shown for HYPO vs. INIT (*), RLD vs. HYPO ($) (see Supplementary Table5for details).The number of symbols indicates the statistical level (three for p < 0.001).

Figure 8 .
Figure 8. Correlation analysis of relative adjustments from INIT to HYPO.(a) Correlations between relative adjustments in spatiotemporal parameters of running and in motor primitives, as measured by linear metrics.(b) Correlations between relative adjustments in spatiotemporal parameters of running and in motor primitives, as measured by nonlinear metrics.Spearman's rank correlations were performed on the data from the participants presenting the synergy of interest (see Supplementary Table3).The blue shaded areas correspond to the 95% confidence interval.Non-significant correlations are shown in Supplementary Table6.

Table 5
for details).The number of symbols indicates the statistical level (three for p < 0.001).Vol.:(0123456789) Scientific Reports | (2024) 14:2707 | https://doi.org/10.1038/s41598-023-50076-6 wrapped with a self-adhesive and a tubular bandage around the thigh and shank.The recorded force and sEMG signals were synchronized.All data analyses were performed on the last 60 running cycles of each condition.